home *** CD-ROM | disk | FTP | other *** search
/ Fritz: All Fritz / All Fritz.zip / All Fritz / FILES / PROGMISC / PCSSP.LZH / PC-SSP.ZIP / MATINSL.ZIP / SINV.FOR < prev   
Text File  |  1985-11-29  |  4KB  |  121 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE SINV
  5. C
  6. C        PURPOSE
  7. C           INVERT A GIVEN SYMMETRIC POSITIVE DEFINITE MATRIX
  8. C
  9. C        USAGE
  10. C           CALL SINV(A,N,EPS,IER)
  11. C
  12. C        DESCRIPTION OF PARAMETERS
  13. C           A      - UPPER TRIANGULAR PART OF THE GIVEN SYMMETRIC
  14. C                    POSITIVE DEFINITE N BY N COEFFICIENT MATRIX.
  15. C                    ON RETURN A CONTAINS THE RESULTANT UPPER
  16. C                    TRIANGULAR MATRIX.
  17. C           N      - THE NUMBER OF ROWS (COLUMNS) IN GIVEN MATRIX.
  18. C           EPS    - AN INPUT CONSTANT WHICH IS USED AS RELATIVE
  19. C                    TOLERANCE FOR TEST ON LOSS OF SIGNIFICANCE.
  20. C           IER    - RESULTING ERROR PARAMETER CODED AS FOLLOWS
  21. C                    IER=0  - NO ERROR
  22. C                    IER=-1 - NO RESULT BECAUSE OF WRONG INPUT PARAME-
  23. C                             TER N OR BECAUSE SOME RADICAND IS NON-
  24. C                             POSITIVE (MATRIX A IS NOT POSITIVE
  25. C                             DEFINITE, POSSIBLY DUE TO LOSS OF SIGNI-
  26. C                             FICANCE)
  27. C                    IER=K  - WARNING WHICH INDICATES LOSS OF SIGNIFI-
  28. C                             CANCE. THE RADICAND FORMED AT FACTORIZA-
  29. C                             TION STEP K+1 WAS STILL POSITIVE BUT NO
  30. C                             LONGER GREATER THAN ABS(EPS*A(K+1,K+1)).
  31. C
  32. C        REMARKS
  33. C           THE UPPER TRIANGULAR PART OF GIVEN MATRIX IS ASSUMED TO BE
  34. C           STORED COLUMNWISE IN N*(N+1)/2 SUCCESSIVE STORAGE LOCATIONS.
  35. C           IN THE SAME STORAGE LOCATIONS THE RESULTING UPPER TRIANGU-
  36. C           LAR MATRIX IS STORED COLUMNWISE TOO.
  37. C           THE PROCEDURE GIVES RESULTS IF N IS GREATER THAN 0 AND ALL
  38. C           CALCULATED RADICANDS ARE POSITIVE.
  39. C
  40. C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  41. C           MFSD
  42. C
  43. C        METHOD
  44. C           SOLUTION IS DONE USING THE FACTORIZATION BY SUBROUTINE MFSD.
  45. C
  46. C     ..................................................................
  47. C
  48.       SUBROUTINE SINV(A,N,EPS,IER)
  49. C
  50. C
  51.       DIMENSION A(1)
  52.       DOUBLE PRECISION DIN,WORK
  53. C
  54. C        FACTORIZE GIVEN MATRIX BY MEANS OF SUBROUTINE MFSD
  55. C        A = TRANSPOSE(T) * T
  56.       CALL MFSD(A,N,EPS,IER)
  57.       IF(IER) 9,1,1
  58. C
  59. C        INVERT UPPER TRIANGULAR MATRIX T
  60. C        PREPARE INVERSION-LOOP
  61.     1 IPIV=N*(N+1)/2
  62.       IND=IPIV
  63. C
  64. C        INITIALIZE INVERSION-LOOP
  65.       DO 6 I=1,N
  66.       DIN=1.D0/DBLE(A(IPIV))
  67.       A(IPIV)=DIN
  68.       MIN=N
  69.       KEND=I-1
  70.       LANF=N-KEND
  71.       IF(KEND) 5,5,2
  72.     2 J=IND
  73. C
  74. C        INITIALIZE ROW-LOOP
  75.       DO 4 K=1,KEND
  76.       WORK=0.D0
  77.       MIN=MIN-1
  78.       LHOR=IPIV
  79.       LVER=J
  80. C
  81. C        START INNER LOOP
  82.       DO 3 L=LANF,MIN
  83.       LVER=LVER+1
  84.       LHOR=LHOR+L
  85.     3 WORK=WORK+DBLE(A(LVER)*A(LHOR))
  86. C        END OF INNER LOOP
  87. C
  88.       A(J)=-WORK*DIN
  89.     4 J=J-MIN
  90. C        END OF ROW-LOOP
  91. C
  92.     5 IPIV=IPIV-MIN
  93.     6 IND=IND-1
  94. C        END OF INVERSION-LOOP
  95. C
  96. C        CALCULATE INVERSE(A) BY MEANS OF INVERSE(T)
  97. C        INVERSE(A) = INVERSE(T) * TRANSPOSE(INVERSE(T))
  98. C        INITIALIZE MULTIPLICATION-LOOP
  99.       DO 8 I=1,N
  100.       IPIV=IPIV+I
  101.       J=IPIV
  102. C
  103. C        INITIALIZE ROW-LOOP
  104.       DO 8 K=I,N
  105.       WORK=0.D0
  106.       LHOR=J
  107. C
  108. C        START INNER LOOP
  109.       DO 7 L=K,N
  110.       LVER=LHOR+K-I
  111.       WORK=WORK+DBLE(A(LHOR)*A(LVER))
  112.     7 LHOR=LHOR+L
  113. C        END OF INNER LOOP
  114. C
  115.       A(J)=WORK
  116.     8 J=J+K
  117. C        END OF ROW- AND MULTIPLICATION-LOOP
  118. C
  119.     9 RETURN
  120.       END
  121.